Measuring orbital eccentricity and periastron advance in quasi-circular black hole 

simulations. 
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We compare different methods of computing the orbital eccentricity of quasi-circular binary black 
hole systems using the orbital variables and gravitational wave phase and frequency. For eccen- 
tricities of about a per cent, most methods work satisfactorily. For small eccentricity, however, the 
gravitational wave phase allows a particularly clean and reliable measurement of the eccentricity. 
Furthermore, we measure the decay of the orbital eccentricity during the inspiral and find reasonable 
agreement with post-Newtonian results. Finally, we measure the periastron advance of non-spinning 
binary black holes, and we compare them to post-Newtonian approximations. With the low uncer- 
tainty in the measurement of the periastron advance, we positively detect deviations between fully 
numerical simulations and post-Newtonian calculations. 

PACS numbers: 04.25. D-, 04.25.dg, 04.25.Nx, 04.30.-w, 04.30.Db 



I. INTRODUCTION 

The inspiral and merger of binary black holes or neu- 
tron stars is one of the most promising sources for cur- 
rent and future generations of gravitational wave detec- 
tors such as LIGO and VIRGO. The late stage of the 
inspiral, corresponding to the final few orbits and merger 
of the binary, is highly dynamical and involves strong 
gravitational fields, and it must be handled by numeri- 
cal relativity. Breakthroughs in numerical relativity have 
allowed a system of two inspiraling black holes to be 
evolved through merger and the ringdown of the rem- 
nant black hole pHl5| . 

During the inspiral of an isolated binary, the orbit cir- 
cularizes via the emission of gravitational waves flfjl . 
As a result, even binaries starting with some eccentricity 
at the beginning of their stellar evolution are expected to 
have negligible eccentricity by the time the frequency of 
the emitted gravitational radiation enters the frequency 
band of ground based detectors. 



However, different physical scenarios [18H26I ] suggest 
that binaries could approach merger with a significant 
eccentricity without being circularized by radiation reac- 
tion. This implies that eccentric binaries are a poten- 
tial gravitational wave source for ground based interfer- 
ometers. For example, in globular clusters, the Kozai 
mechanism (l8j could increase the eccentricity of an in- 
ner binary's orbit through a secular resonance caused 
by a third perturbing black hole on an outer orbit [20| . 
Many-body encounters of black holes in globular clusters 
could also result in the merger of highly eccentric bina- 
ries Ref. 24 1 predicted that 30% of the hierarchical 



triple black hole systems formed in a globular cluster will 
possess eccentricities greater than 0.1 when their emitted 
gravitational waves pass through a frequency of 10Hz. 

For these reasons considerable attention has been paid 
to eccentric binaries. Analytical waveform templates 



have been constructed for the gravitational wave signal 
emitted by compact binaries moving in inspiraling eccen- 
tric orbits |27H29| . In this case, orbits involve three dif- 
ferent time scales: orbital period, periastron advance and 
radiation reaction time scales. By combining these three 
time scales, one computes "postadiabatic" short-period 
contributions to the orbital phasing and gravitational 
wave polarizations. These gravitational wave polariza- 
tions are needed for astrophysical measurements with 
gravitational wave interferometers. Refs. j3fj| - l32l | inves- 
tigated the impact of eccentricity on gravitational wave 
detection, specifically the potential loss in the signal-to- 
noise ratio when "circular" waveform templates are ap- 
plied to search for eccentric binaries. 

Eccentric black hole binaries have also been studied 
with direct numerical simulations. Ref. [33| studied the 
variation of the signal to noise of the eccentric evolutions 
of intermediate mass binary black hole mergers as a func- 
tion of mass and eccentricity. Ref. }34| presented binary 
black holes in zoom-whirl orbits where the waveforms 
are modulated by the harmonics of these zoom-whirls. 
In Ref. |35[, the authors studied the transition from in- 
spiral to plunge in general relativity by computing grav- 
itational waveforms of eccentric nonspinning, equal mass 
black-hole binaries. They analyzed the radiation of en- 
ergy and angular momentum in gravitational waves, the 
contribution of different multipolar components and the 
final spin of the remnant black hole. Ref. [HI presented 
results from numerical simulations of equal-mass, non- 
spinning binary black hole inspiral and merger for various 
eccentricities, and they measured the final mass and spin 
of the remnant black hole. Ref. 37J compared a numer- 
ical relativity simulation of an eccentric binary system 
with eccentricity 0.1 with corresponding post-Newtonian 
(PN) results. They found better agreement when the 
eccentric PN expressions are expanded in terms of the 
frequency-related parameter x = (flM) 2 ^ 3 , where f2 is or- 
bital frequency and M is total mass of the binary, rather 
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than the mean motion n = 2tt/ P, where P is the orbital 
period. 



Beyond the Newtonian limit, the orbital eccentricity 
is not uniquely defined and a variety of definitions have 
appeared in the literature. Ref. [H[ used a definition of 
the eccentricity for which a Newtonian orbit is momen- 
tarily tangent to the true orbit (the "osculating" eccen- 
tricity) , while other authors [27H29L l39j defined multiple 
"eccentricities" to encapsulate different aspects of non- 
circular orbits at PN order. Another useful definition 
for large eccentricity in numerical simulations is given in 
Refs. ESEI. 



Similarly, numerical relativists [42H46J introduced sev- 
eral methods for defining and measuring the eccentricity 
using the residual oscillations in the orbital frequency, 
proper horizon separation and coordinate separation. 
These eccentricity definitions are necessary to compare 
the numerical waveforms with the waveforms produced 
by analytic techniques (i.e., PN methods). They behave 
differently depending on the magnitude of the eccentric- 
ity and details of the numerical simulation, like employed 
gauge conditions, or presence of numerical noise. This 
makes it important to specify the validity regimes of these 
definitions. 



This paper deals with two related topics: First, we re- 
visit many of the eccentricity definitions used so far in 
numerical work and compare them systematically. Wc 
find that for eccentricities of a few percent, most defi- 
nitions work satisfactorily. However, for very small ec- 
centricities, e ~ 10~ 4 , computation of the eccentricity 
based on the extracted gravitational waves is superior. 
In the second part of the paper, we measure decay of or- 
bital eccentricity and periastron advance for inspiraling 
black hole binaries, and compare these measurements to 
post-Newtonian calculations. 



Section HT1 summarizes eccentricity definitions that are 
useful for measuring eccentricity in quasi-circular runs. 
In Section [TTl we compare these approaches, as well as 
some new ones, on the 15-orbit inspiral presented by 
Boyle et al. 14| and on the data of a new simulation 
of an eccentric (e = 0.05) nonspinning equal mass binary 



black hole. Next, by measuring the extrema in the eccen- 
tricity estimator, we estimate in Sec. IHII the decay of the 
eccentricity of these runs as well as the radial frequency. 
This allows us in Section IIVI to estimate the periastron 
advance for these runs from the ratio of the orbital fre- 
quency to the radial frequency as well as the periastron 
advance of a set of quasi-circular nonspinning binaries of 
mass ratios 2, 3, 4 and 6. The numerically estimated pe- 
riastron advance is then compared to the 3PN formula of 
the periastron advance [2^, [5^, [47| ■ 



II. ECCENTRICITY ESTIMATORS 

A. Definitions 

For a non-precessing binary in an orbit with zero eccen- 
tricity, orbital variables and their time derivatives change 
monotonically as the holes inspiral to merger. In numer- 
ical simulations, however, a small eccentricity is intro- 
duced by imperfections of the initial data. As a result, 
small residual oscillations with amplitude proportional to 
the eccentricity are added to the monotonically changing 
orbital variables and their derivatives. To estimate the 
eccentricity, one needs to determine these residual oscil- 
lations. 

Different methods to estimate the eccentricity [42], [iil . 
|45| used the orbital frequency, separation between the 
holes (coordinate or proper separation), or some New- 
tonian formula containing both of these variables. Sim- 
ilarly, time derivatives of these variables could be used 
in these definitions of the eccentricity. Basically all ap- 
proaches construct an eccentricity estimator exit) such 
that for Newtonian orbits 



ex(t) = e cos(fl r t + 



(1) 



where e is the eccentricity 1 and f2 r is the frequency of 
radial oscillations in the quasi-circular orbit. The key 
property of ex (t) is that it is an oscillating function with 
amplitude equal to e. 

In order to define eccentricity for general relativistic in- 
spirals, one computes a tentative eccentricity estimator 
ex(t), and checks its behavior. If it behaves as Eq. ((TJ), 
one reads off the eccentricity e as the amplitude of the 
oscillations. The resulting eccentricity estimates are not 
local in time nor continuous functions of time but rather 
orbit-averaged quantities. Deviation from sinusoidal be- 
havior indicates that particular eccentricity estimator is 
not reliable, and one must verify to what extent the ec- 
centricity estimators behave as expected and to what ex- 
tent they agree. 

The estimated value of the eccentricity will differ 
slightly depending on the method used and the noise in 
the numerical data. In this paper, we compare typical 
eccentricity estimates using a Newtonian formula as in 
Ref. 1421 or the orbital frequency and separation as in 
Ref. [45j . These eccentricities are also compared to new 
ones computed from the wave phase and frequency ex- 
tracted at a given radius. Other definitions of the eccen- 
tricity could be used, but we restrict the study to these 
typical definitions. 

To make this rather abstract discussion more concrete, 
consider the Newtonian formula for the radial distance d 
between the two objects with eccentricity eNcwt 

d(t) = d [1 + e Nc wt cos(ft r t + O )] + 0{e 2 ) . (2) 



1 The eccentricity e is well-defined for Newtonian orbits. 
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Based on this formula, one can define the eccentricity 
estimator e^i) 



d(t) 



d(t) 



ecos(fi r i + </>o), 



(3) 



where the average distance d equals do in Newtonian 
gravity. For a general relativistic system, one obtains 
d(t) by a fit over several radial oscillation periods. If the 
residual d(t)—d(t) oscillates sinusoidally — which it indeed 
does for sufficiently large eccentricity — the amplitude of 
these oscillations defines an associated eccentricity ej- 

From the trajectory of the two objects, one can also use 
the orbital phase and frequency to define the correspond- 
ing eccentricity estimators using the following Newtonian 
relation 0: 
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<f> = M + 2esmM + -e 2 sm2M + 0{e 3 ), (4) 

where M. is the mean anomaly and <f> is the orbital phase. 
Equivalent to Eq. Q for numerical simulations is the 
relationship 



$(i) = $ + tt Q t + 2esin(^) + 0(e 2 ) . 



(5) 



where O is the average fitted orbital frequency and $o is 
some phase offset. Then the eccentricity estimator e$(t) 
is written as 



5*(t) = 



(6) 



Prom the time derivative of Eq. ([FJ) and the replacement 
we obtain an eccentricity estimator in terms 
of the orbital frequency (as in Ref. |45j ) 



en(t) 



n(t) - n 



(7) 



Notice that since the radial oscillation results from ec- 
centricity, fi r is different from f2 , the average of the or- 
bital frequency. The eccentricities of Eqs. © and J7]) 
will differ by a factor f2 r /Oo- For Newtonian orbits, 
f2 r /f2o = lj and this factor drops out. But for the bi- 
nary black hole case, the factor is about 1.4, causing the 
difference between Figs. [4] and [5] below. This is easily 
seen by writing the eccentricity estimator from Eqs. ([5]) 
and (JT]) as 



2n 



= eQr/Qo sin(f2 r t) . 



(8) 



FIG. 1: The equal mass nonspinning binary run with eccen- 
tricity e ~ 0.05. As a function of time, the top panel shows 
the proper horizon separation and the bottom panel shows 
the orbital frequency. For such a value of the eccentricity, it 
is easy to measure the decay rate of the eccentricity and esti- 
mate the periastron advance of the binary near the merger. 



We will primarily analyze the 16 orbit long inspiral sim- 
ulation of an equal mass, non-spinning black hole binary 
presented in Ref. |l4| (specifically, the run labeled 30c-l). 
This run with eccentricity of about 6 x 10~ 5 is used to 
compute the eccentricity data in Figs. El [3j H] and El To 
compute eccentricity estimators, we use the orbital fre- 
quency f2, the coordinate separation between the holes 
D, the proper horizon separation s (defined as the inte- 
grated distance between the holes along the coordinate 
axis, cf. Ref. [lH) as well as the gravitational wave phase 
4> and the gravitational wave frequency to. 

Furthermore, we utilize recent runs of quasi-circular 
nonspinning binaries [50( | with mass ratios 2, 3, 4 (last- 
ing 15 orbits) and mass ratio 6 (lasting 8 orbits). The 
eccentricity of these runs is also of the order of mag- 
nitude 10 -5 . The periastron advance and the resulting 
frequency modulation are estimated in Fig. [51 

As a separate check, another equal mass nonspinning 
binary with moderate eccentricity (e ~ 0.05) is evolved to 
compare various eccentricity estimators and measure the 
periastron advance for a case that is not quasi-circular. 
Figure [T] shows the proper separation as well as the or- 
bital frequency as a function of time for this eccentric 
binary. 



C. A Newtonian definition 



B. Numerical data 

Before introducing several further eccentricity estima- 
tors, let us briefly describe the numerical binary black 
hole simulations that we will analyze. All runs have been 
performed with the Spectral Einstein Code (SpEC) (49| . 



The first use of eccentricity estimators was by Buon- 
nano, Cook & Prctorius who consider the following 
relationship that holds for Newtonian orbits with eccen- 
tricity e No wt: 

[^{tfritf/M - 1] = e Nowt cos <j>{t) . (9) 
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FIG. 2: Eccentricity estimator eecp [42| applied to a simu- 
lation with e ~ 0.05 (top panel) and e ~ 6 x 10~ 5 (bottom 
panel). The dashed and solid lines correspond to eBCp(i) com- 
puted from the coordinate separation and the proper horizon 
separation. For the large eccentricity run, escp exhibits clear 
oscillations, whereas for the small eccentricity run, eBCP is 
dominated by other features. In both cases, the amplitude of 
eBCP is smaller when defined using coordinate distance D. 



Here fl^,(t) and <j)(t) denote orbital frequency and phase, 
respectively, and r is the separation of the masses. Mo- 
tivated by Eq. ([9]), Buonnano, Cook & Prctorius define 
an eccentricity estimator 

e BC p(t) = il+(t) 2 r(t) 3 /M - [n^t) 2 r(t) 3 /M] Rt , (10) 

where now fi^i) and r(t) are extracted from the numer- 
ical simulation. To compute this eccentricity estimator 
eBCPi we fit the function ^^(tfrft) 3 /M to a polynomial 
in time, 



(11) 



i=0 



We found that a fifth order polynomial ensures a good fit. 
The polynomial order needs to be high enough to reliably 
capture the smooth inspiral trend in n(t) 2 r(t) 3 /M, but it 
should not capture the higher frequency oscillations due 
to eccentricity. When applying this procedure to a binary 
black hole inspiral, one has to decide how to generalize 
the Newtonian separation r(t) to curved space. We use 
two choices, the coordinate distance D{t) between the 
centers of the apparent horizons, and the proper separa- 
tion s(t) between the apparent horizons, computed along 
a straight coordinate line connecting the centers of the 
apparent horizons. 

In Fig. [3J we plot the eccentricity estimator eecp com- 
puted using the coordinate separation and proper horizon 
separation as described above. In the top panel, we plot 



ep>cp(£) using the binary run with eccentricity e ~ 0.05. 
Using the proper horizon separation s, the estimated ini- 
tial eccentricity, 0.07, is larger by nearly a factor of 2 
than in the case where the coordinate separation D is 
used (0.03). This is due to different numerical values 
for the distances, (s/D) 3 <~ 1.8. Both eccentricity es- 
timators are in phase during the whole time interval as 
expected. In both cases, the eccentricity magnitude de- 
creases between t = and t = 2500M. In this case, a 
clear decaying sinusoidal signal is obtained without any 
higher harmonics showing up at later times. 

In the bottom panel, we examine the equal mass bi- 
nary with eccentricity e~6x 10 -5 . For this case, no 
clean sinusoidal signal is apparent. While ep,CP com- 
puted from s(t) shows oscillations, they are faster than 
the orbital period, and can therefore not be attributed to 
orbital eccentricity. Because eecp does not show the ex- 
pected behavior, it is not meaningful to attribute a value 
of eccentricity to this analysis. For these small eccen- 
tricities, eecp is dominated by other effects, possibly the 
coordinate dependence of the separation measurements. 



D. Eccentricity from orbital variables 

Husa et al [45j fitted directly the orbital frequency fi(i) 
or the coordinate separation D(t) to a function of the 
form 



i/2 



(12) 



with fitting parameters t m , the coalescence time, and the 
coefficients ai. The eccentricity estimator is then defined 
as 



ex(t) 



^nr(£) — Xft t (t) 
kX &t (t) 



(13) 



where Xnr(^) is the numerical orbital variable and Xfa(t) 
is the polynomial fit of -X"nr (£)• We shall compute three 
eccentricity estimators using Eq. (|12p . which differ in the 
quantity being fitted: e s (t) and e d (t) are based on proper 
separation and coordinate separation between the black 
holes, with the value k = 1; en(t) uses the orbital fre- 
quency, where k = 2. In the Newtonian limit, these 
estimators are identical to first order in eccentricity. 

Figure [3] shows these eccentricity estimators for a run 
with fairly large eccentricity and for a run with very small 
eccentricity. For large eccentricity e = 0.05, the various 
eccentricity estimators have a smooth decaying sinusoidal 
signal. This allows measuring a nearly identical value of 
the eccentricity for the three orbital variables from the 
amplitude of the residual oscillations. The phasing is also 
consistent between the different eccentricity estimates: 
The orbital frequency is a maximum when the separation 
is a minimum and vice-versa. 

In the bottom panel of Fig. [31 we plot the eccentricity 
estimators applied to a simulation with much smaller ec- 
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FIG. 3: Eccentricity estimators based on orbital trajecto- 
ries applied to simulations with eccentricity e ~ 0.05 (upper 
panel) and e~6x 10 _J (lower panel). The quantities en, e s 
and eu are computed from orbital frequency, proper horizon 
separation and coordinate separation using Eq. (|13[) . 



ccntricity e~6x 10 -5 . The behavior of ejj and en is er- 
ratic. Higher-order harmonics are clearly visible, and the 
cxtrcma are not monotonically decreasing, as one would 
expect from the circularizing effect of gravitational radi- 
ation. However, e s shows no increase in the eccentricity 
during the late stages of the inspiral, and no additional 
significant harmonics appears even at t = 3500 M. The 
order of the polynomial fit depends on the time range of 
the fit. In this case, a fifth order polynomial was enough 
to capture the oscillatory behavior in the eccentricity es- 
timator in the time range 500 M < t < 3500 M. Note that 
the orbital phase could also be used to measure the ec- 
centricity estimator using Eq. (fl"3"| (but without division 
by X fit ). 



E. Eccentricity from gravitational waves 

All eccentricity estimators discussed so far utilize 
coordinate-dependent quantities like separation or or- 
bital frequency. Therefore, one might suspect that the 
higher harmonics visible in Figs. [5] and are caused by 
gauge effects. The gravitational radiation at future null 
infinity is expected to be gauge-invariant, removing the 
dependence on gauge-dependent quantities. These con- 
siderations motivate the use of the gravitational wave 
phase and frequency to define eccentricity. 

We extract the (I, m) = (2, 2) mode of the gravitational 
wave using the Newman-Penrose scalar ^4 and define the 



wave phase <f>(t) as [ijj 

*f (r,t) = A(r,t)e- i0(r ' 4) . (14) 
Then the gravitational-wave frequency is defined as 

rlrh 

(15) 



dt 



The waveforms extracted at finite radii are extrapolated 
to null infinity using the procedure in j5l| . The wave 
phase <p and frequency uj are measured as a function of the 
retarded time t — r* , where r* is the tortoise-coordinate 
radius defined as 



2M A dm In 



\2M 



- 1 



DM 



(16) 



where Madm is the ADM mass of the initial data. At 
early times, the gravitational waveforms are contami- 
nated by high frequency noise from imperfect initial data. 
To measure the amplitudes and locations of the extrema 
in the eccentricity estimator more accurately, the residual 
functions are filtered using a low-pass Butterworth filter 
with the Matlab function "filtfilt" [Hj]. The filtered data 
can be used to measure the eccentricity for retarded time 
t-r*> 1000M. 

Based on the gravitational wave phase, we define the 
eccentricity estimator 



>(*) = 



4>NB.(t) - 0fit(*) 



(17) 



where an additional factor of 1/2 arises because the wave 
phase is approximately twice the orbital phase. 

In Fig. |H we plot the eccentricity estimator computed 
from the gravitational wave phase of the (2,2) mode ex- 
tracted at the radii r = 75M, r = 240M and extrapo- 
lated to infinity using terms up to 1 /r 2 versus t — r*. The 
eccentricity estimate is independent of the radius value at 
which the wave is extracted, and various estimates agree 
to within 5% in both amplitude and phase for different 
radii of extraction. 

Using the wave frequency we define the eccentricity 
estimator e^it) 



,{t) 



2wat(*) 



(18) 



Computation of the gravitational wave frequency u = 
d<p/dt requires a derivative of <f>(t), which increases nu- 
merical noise. Given the small amplitude of the ef- 
fect under consideration (the fractional change in uj is 
2e = O(10~ 4 )), the increased noise noticably affects e w . 
It is usable only at finite extraction radius, and even there 
only for t - r* > 2000M. 

In Fig. [SJ we compute the eccentricity estimator from 
the wave frequency extracted at r = 75M and r = 240 M. 
The extrapolated data to infinity is not shown because of 
its sensitivity to noise. The two curves have a nearly si- 
nusoidal behavior with the phase agreeing to within 10%. 
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FIG. 4: Eccentricity estimator computed from the gravita- 
tional wave phase as a function of retarded time t — r*. In this 
plot, the eccentricity estimator is computed from the gravita- 
tional wave extracted at finite radii r — 75M and r — 240A/ 
and from data extrapolated to infinity. The three curves agree 
in amplitude and phase to within 5% in the retarded time in- 
terval 1000M < t - r* < 3000Af. 



FIG. 5: Eccentricity estimator e u computed from the grav- 
itational wave frequency as a function of the retarded time 
t — r*. In this plot, the eccentricity estimator is computed 
from the gravitational wave extracted at r — 75M and 
r — 240M. The eccentricity estimator is contaminated by 
significant noise caused by imperfect initial data at a time 
earlier than t/M = 2000. 



However, the amplitude differs by 25% between the wave 
data measured at r = 75M and r = 240M. The reduced 
sensitivity to noise is an important advantage of over 

For the binary with eccentricity 0.05, plots similar to 
Figs. @] and [3] with smooth sinusoidal behavior could eas- 
ily be obtained. 

Computation of the eccentricity from gravitational ra- 
diation (e^ and ej) is better behaved than the methods 
using orbital variables. Only one harmonic mode ap- 
pears in the data — even for the low-eccentricity run with 
e~6x 10~ 5 — and the eccentricity is decreasing as the 
binaries inspiral toward merger. We attribute this im- 
provement to the disappearance of coordinate and gauge 
effects when the data are extracted further away from 
the holes. 

The eccentricities extracted from and e w in Figs. @] 
and [5] are inconsistent with each other; they differ by a 
factor fi-r/fi^ as explained in Sec. Ill Al 

One might also consider a definition of the eccentricity 
based on taking the time derivative of the wave frequency. 
From Eq. the second time derivative of the orbital 
phase is given by: 

<i> = M -2e{M cos M+M 2 sin M) + 0{e 2 ), (19) 

where the amplitude of the oscillatory part is 

2ev A4 2 + M. A . Then, the eccentricity estimator com- 
puted from the time derivative of the wave frequency e<f w 



is then defined as 



<?W - 0fit 



4 



(20) 



The main advantage of a such a definition is that it re- 
quires a lower order fitting polynomial. Unfortunately, 
the numerical derivatives necessary to compute <j) am- 
plify noise, and so this method becomes impractical for 
the numerical evolutions considered. 

In Table HI we summarize the eccentricity definitions 
examined in this paper, the data range between U/M and 
tf/M employed in the fits, and the order of the fitting 
polynomial n for the 15-orbits quasi-circular nonspinning 
binary. We also give an estimate of the eccentricity value 
at t/M = 1000,2000 and 3000 and its estimated error 
Se/e for each method. 



III. BEHAVIOR OF ECCENTRICITY DURING 
INSPIRAL 

Radiation reaction reduces eccentricity during the in- 
spiral of a binary compact object, as shown by the 
post-Newtonian calculation by Peters (Itij . Using the 
quadrupole approximation, Peters derived the evolution 
of the orbital eccentricity during the inspiral caused by 
the emission of gravitational waves. In the limit of small 
eccentricity, the eccentricity is related to the semi-major 
axis a by 



e oc a 



19/12 



(21) 
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Aw/(2oj flt ) 


1922 


3861 


7 


- 


4.3 KlO -6 


3.7xl0^ 5 


15-25% 


Coordinate distance 


eo 


AI)/J) fi t 


480 


3367 


7 


6.7xl0~ 5 


4.9xl0~ 5 


6.3xl0~ 5 


15-40% 


Proper separation 


eh 




480 


3367 


5 


5.0 xl0~ 5 


3.9 xl0 _s 


3.4xl0 -5 


10-20% 


Orbital frequency 


en 


Afi/(2fi fit ) 


480 


3367 


7 


6.2xl0 -5 


4.1xl0 -5 


3.4xl0 -5 


20-30% 


BCP 


eBCP 


A(fi(i) 2 r 3 ) 


480 


3367 


5 


3.5xl0 -5 


2.4xl0" 5 


2 xlO" 5 


50-80% 



TABLE I: Summary of the eccentricity measurement methods. U (tf) is the initial (final) time of fitting, n is the order of the 
best fitting polynomial in the time interval [ti/M,tf/M\. e is the eccentricity estimate at the time t with the relative error 
Se/e. 



CD 
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0.07 
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0.04 - 

0.03 
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□ 1000c. (low ecc. nan) 
O lOOOe (low ecc. nan) 

(high ecc. run) 
x e (high ecc. nan) 



10 11 12 13 14 15 16 

s/M 

FIG. 6: Eccentricity as a function of proper horizon separa- 
tion. We show data for two simulations, with high and low 
eccentricity. For each run we compute eccentricity from the 
GW-phase <j> an d the proper separation s. The dashed line 
represents the power-law s 19 ^ 12 predicted by post-Newtonian 
theory (See. Eq. l2Tj) . 



to maximum, or vice versa) of the eccentricity estimator 

e = |Amin 2 Amax| . (22) 

We further associate this eccentricity with the time half- 
way between the two extrema under consideration: 



t(e) 



t(A m i n ) + i(A max ) 



(23) 



At the time of this average eccentricity, the separation 
is measured numerically. In the case when gravitational 
wave data is used, the wave phase is approximated as 
a function of the separation by using the retarded time 
t - r*. 

The results are plotted in Fig. |6] Fitting a power-law 



log e = a + p log s 



(24) 



to the numerical data yields j3 w 1.4. These decay esti- 
mates are in reasonable agreement with Peters' predic- 
tion (/3 = 19/12 ~ 1.583), as can be seen by the indicated 
power-law in Fig. [6] The orbital eccentricity decays sim- 
ilarly in the two simulations with different eccentricity. 



The first confirmation of the decay of eccentricity in a 
fully numerical binary black hole inspiral was presented 
by Pfeiffer et. al. [44|. Pfeiffer et. al. measured the decay 
rate of the eccentricity for an equal mass, nonspinning 
binary with an eccentricity of about 0.02 during the last 
five orbits of the inspiral. The precise decay rate de- 
pended on the definition of the eccentricity used. For a 
definition based on the orbital frequency, good agreement 
with Eq. (f2~T|) was found. 

In Sec. mi we established that the eccentricity estima- 
tors (wave phase) and e s (proper horizon separation) 
show the cleanest oscillatory behavior. Using these two 
eccentricity estimators, we compute as follows the eccen- 
tricity as a function of time for the much longer inspirals 
considered here. We first define the "average" eccen- 
tricity over one half of a radial oscillation as the differ- 
ence between two consecutive extrema (from minimum 



IV. PERIASTRON ADVANCE 

The pcriastron advance is one of the new features for 
rclativistic eccentric orbits that is not present in Newto- 
nian gravity. It has been computed analytically in the 
post-Newtonian regime up to third order but — to our 
knowledge — it has never been estimated numerically in 
binary black hole simulations. Periastron advance will 
lead to a modulation of the gravitational wave signal for 
eccentric binaries and will impact gravitational wave de- 
tection strategies. Therefore, it is important to know 
what this frequency is and how it changes as a function 
of the mass ratio. The fractional periastron advance per 
orbit, K, is defined as 



where A$ = $ — 2tt is the periastron advance per orbit. 
The dimensionless parameter K is related to the radial 
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frequency f2 r and the orbital frequency J7$ through 

1 . (26) 



1.6 



Numerical method for measuring the periastron 
advance 



From an eccentricity estimator ex, cf. Eq. (|13|) . one 
can read off not only the eccentricity (via the amplitude 
of ex), but also the frequency of the radial motion, tt r 
(from the oscillation period). We shall define the period 
of the radial oscillation as twice the time interval between 
two consecutive extrema (from minimum to maximum, or 
vice versa) in the eccentricity estimator curve. We em- 
ploy the following procedure to compute the periastron 
advance: 

f . Choose a cleanly oscillating eccentricity estimator 
ex(t). We will use e^, cf. Fig. 21 

2. Find the extrema of ex{t)- This gives a time list 
(to, t\, tk, • ••) corresponding to all perihelia or 
aphelia (i.e., extrema in the residual radial veloc- 

ity). 

3. Interpolate the orbital phase $ to the times tk- Be- 
tween neighboring data points, the orbital phase 
changes by $(tk+i) — $(tfc_i), whereas the ra- 
dial phase changes by 2ir. Therefore, the ra- 
tio between orbital and radial phase increase is 
($(i fe+ i) - $(t fe _i))/27T, and so 



fi<j> 



$(t 



k+l, 



*(*k-l) 



2- 



(27) 



For the very low eccentricity simulation (e ~ 5 x 10 -5 ), 
the periastron advance is very difficult to measure be- 
cause the amplitude of is so small. The uncertainty in 
the extracted Q,^,/Vl r is about ±0.1 for 0.02 < Mf2$ < 
0.03. The error in the estimated periastron advance in- 
creases at higher frequencies as the binary evolves closer 
to merger. The eccentricity estimators depend on details 
of the polynomial hts, and it is more difficult to read off 
these small eccentricity residuals near the plunge. There- 
fore, f2$/f2 r for the e ~ 5 x 10~ 5 run should not be 
trusted for Mfl^ > 0.03. 

In the simulation with larger eccentricity e ~ 0.05, 
by contrast, the periastron advance is easier to measure, 
because the amplitude of e$ is proportional to the eccen- 
tricity. We obtain correspondingly smaller errors, about 
3% at frequencies Mf2 < 0.03. While we are able to ex- 
tract 57$ /VL r at higher frequencies for the simulation with 
e ~ 0.05, recall that the numerical data are constructed 
from consecutive extrema of e^. At late times (close to 
merger), there is an increasing amount of orbital evo- 
lution during such an interval, which renders ambiguous 
both the definition of / Q r and its association with one 
orbital frequency. 



a" 



1.2 



o 


q=l, e~0.05 


X 


q=l, e~5xl0" 5 




Test-mass 




q=l, 3PN 




0.02 



0.025 

MQ. 



0.03 



<f> 



FIG. 7: Periastron advance for equal-mass binaries. Plotted 
is the ratio of orbital frequency to radial frequency, Q<s>/Q r , 
versus the orbital frequency A/fi$. The data represent nu- 
merical simulations of equal mass nonspinning black-hole bi- 
naries with two different eccentricities e. Also shown are the 
prediction of post-Newtonian theory for q = 1 and the test- 
mass result based on geodesic motion in Schwarzschild (both 
in the limit e<l). For e~5x 10 -5 , the numerical data is 
unreliable for > 0.03 (see text). 



Figure [7] shows the computed periastron advance for 
the two equal-mass simulations considered here. To facil- 
itate comparison with analytical estimates (see next sec- 
tion), we plot f2$/f2 r as a function of orbital frequency. 
The latter is approximated as half the gravitational wave 
frequency. (This is justified because the deviation from 
this value is much smaller than the error in estimating 
the eccentricity and the periastron advance.) We will 
discuss this figure in the next subsection. 



B. Results 

From Fig.[7]we see that $}$/Q r is positive (i.e. the fully 
general relativistic calculation produces indeed a perias- 
tron advance), and the periastron advance increases with 
increasing orbital frequency MJ7$ , again consistent with 
expectations. The solid and the dashed lines in Fig. [7] 
indicate the periastron advance for a test-mass orbiting 
a Schwarzschild black hole, and for an equal-mass bi- 
nary at 3rd post-Newtonian order (sec Appendix for de- 
tails), and we can now compare these calculations with 
the fully relativistic BBH simulations. The scatter in 
the numerical data fi$/f2 r represents a measure of the 
uncertainty in the periastron advance of the numerical 
simulations. For the e ~ 0.05 simulation, this scatter is 
much smaller than the difference from the 3PN calcula- 
tion. Therefore, we have positively detected a difference 
between fully numeric simulations and 3PN calculations. 
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FIG. 8: Periastron advance for unequal mass BBH. Shown 
is the ratio of orbital frequency to radial frequency, £7$/51 r , 
versus the orbital frequency A/fi<j. for different mass-ratios 
q = Mi/M 2 . 



(f2$/f2 r from the e ~ 10~ 5 simulation coincides with the 
data for the e ~ 0.05 run, although with larger scatter, 
because of trying to extract much smaller variations in 
the numerical data.) The difference between the numeri- 
cal periastron advance curve and the 3PN result is about 
3% at f2$ = 0.02 and continually increases to about 5% 
at 0$ = 0.03. The fully NR periastron advance seems 
to follow more closely the test-mass calculation than the 
equal-mass 3-PN prediction. Note that comparing cither 
of the two analytic results is imperfect: The 3-PN cal- 
culation is for equal masses, but because of the nature 
of post-Newtonian perturbation theory becomes increas- 
ingly less reliable for increasing frequency Mf2$. The 
test-mass limit, in contrast, is an exact calculation, but 
for a system different from an equal-mass binary. Un- 
equal mass binaries with mass ratios very different from 
unity should result in better agreement with the test- 
mass limit, and we will explore this case next. 

Extracting the periastron advance from a series of non- 
spinning unequal mass simulations [5fjj |. we obtain the 
data plotted in Fig. [8] These simulations have very low 
eccentricity in order to accurately model circularized bi- 
naries for gravitational wave data-analysis, with eccen- 
tricities indicated in Fig. [5] The smallness of the eccen- 
tricity is unfortunate for our purposes, as this increases 
the errors in the extracted periastron advance. The pe- 
riastron advance for q = 2 is very similar to the equal 
mass periastron advance data. For higher mass ratio, 
the numerically computed f2$/f2 r seems to increase and 
approach the test-mass result; however, the large uncer- 
tainty in f2$/f2 r for these runs prevents us from drawing 
strong conclusions. 



C. Laplace- Runge-Lenz vector 

The Laplace-Runge-Lenz vector points towards the pe- 
riapsis of the orbit from the center of motion, and there- 
fore it would seem that observing this vector during a 
simulation would result in an immediate measure of the 
periastron advance. This vector is defined in ADM coor- 
dinates in terms of the canonically conjugate position R 
and momentum P as 13911: 



A = P x L- GMfx 



R 
R 



(28) 



where L = R x P and /i is the reduced mass. Unfor- 
tunately, the magnitude of this vector is proportional to 
e, i.e., it will typically be very small. Moreover, it is 
computed as the difference between two large terms that 
almost cancel each other, resulting in large numerical er- 
rors. Furthermore, relativistic effects, such as gauge ef- 
fects, might affect the two terms in Eq. (|28|) differently, 
thus disproportionately affecting the small difference A. 
Yet another obstacle is that the numerical data do not 
give the canonical position and momentum. For all these 
reasons, we found it impossible to measure the perias- 
tron advance from the Laplace-Runge-Lenz vector even 
for the binary run with e <~ 0.05. 



V. DISCUSSION 

We have dealt with three aspects of eccentricity in bi- 
nary black hole simulations: how to measure eccentricity, 
its decay during the inspiral, and periastron advance. 

With regard to techniques to measure eccentricity, this 
paper provides a systematic comparison between several 
different estimators. The ones shown in Figs. [2J [31 H] and 
[S]each displayed a different behavior, even though these 
definitions reduce precisely to the usual eccentricity e in 
the Newtonian limit. Differences appear mainly because 
the data corresponds to a binary in the last phase of 
the inspiral before merger when relativistic effects are 
significant — a regime in which the Newtonian relations 
between the orbital variables are no longer valid. 

The eccentricity estimator eecp (see Fig. [5]) exhibits 
two very undesirable features: For the e ~ 0.05 simula- 
tion, eecp depends strongly on the choice of how sepa- 
ration between the black holes is measured (coordinate 
distance D vs. proper separation s). For small eccen- 
tricities e~5x 10 -5 , no regular oscillatory behavior is 
apparent, rendering escp useless as an eccentricity esti- 
mator. This might be because it uses a definition where 
eccentricity comes in the next-to-leading term, and the 
leading order Newtonian expression is not satisfied. Also, 
the high power of the contribution of the orbital variable 
makes the eccentricity easily affected by high-order har- 
monic modes in the orbital variables. We have observed 
similar behavior when we explored alternative definitions 
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of the eccentricity based on Newtonian formulas combin- 
ing orbital variables. 

Eccentricity measures based on orbital quantities (see 
Fig. [3j give the right amplitude (for t < 2500M in the 
case of en and eo ) , and the phasing is quite consistent be- 
tween the different eccentricity estimators. For instance, 
the orbital frequency is maximal when the separation is 
minimal. However, for the low-eccentricity simulation 
(e~5x 10~ 5 ) higher-order harmonics are clearly visible 
as the binary approaches the merger, in particular for 
the coordinate separation and the orbital frequency 
eo,. The eccentricity measured from the proper horizon 
separation^) is affected least by these coordinate ef- 
fects. 

Eccentricity measures based on extracted gravitational 
waves (see Figs. [4] and [5]) result in clean oscillatory be- 
havior, even for eccentricities as small as considered here. 
No high-order harmonics are noticeable in the wave ex- 
trapolated to infinity during the time interval considered. 
The eccentricity is calculated from the maximum and 
minimum values in the oscillating function without con- 
cern for the coordinate location in the orbit. It is es- 
pecially straightforward to calculate numerically the ec- 
centricity from the wave phase extrapolated to infinity 
without resorting to any notions of "distance" between 
the holes. Computing eccentricity from the gravitational 
wave phase is therefore the preferred method. Unfortu- 
nately, the gravitational wave phase is not as easily acces- 
sible as orbital quantities: One needs to extract gravita- 
tional waves, the waveform is delayed by the light-travel 
time to the extraction radius, and, for best results, one 
may have to extrapolate to infinity. Therefore, in prac- 
tice, eccentricity estimators based on orbital quantities 
may be useful for immediate diagnostics during a simu- 
lation, then confirmed and refined subsequently by eccen- 
tricity estimators based on gravitational wave properties. 

Notice that the eccentricity measurement could be af- 
fected by noise sources such as the "junk radiation" early 
in the simulation or by poor boundary conditions caus- 
ing radiation reflection at the outer boundary. These 
additional oscillations could easily be interpreted as ec- 
centricity. In principle, however, one should be able to 
distinguish them from the eccentricity by the frequency 
of the oscillation. 

The second part of this paper describes measuring the 
decay of orbital eccentricity during the inspiral of equal 
mass non-spinning black hole binaries, revisiting earlier 
work [iij . For both simulations considered, we find that 
eccentricity measured via proper separation (e s ) and via 
gravitational wave frequency decays with the same power 
of proper of separation, s 13 , with exponent /3 « 1.4. This 
is somewhat smaller than the value predicted by post- 
Newtonian expansions, 19/12 1.58. The earlier work, 
which was based on fewer data-points at closer separa- 
tion, found a distinctively smaller exponent when com- 
puting eccentricity from proper separation rather than 
from the orbital frequency. 

The third part of this paper presents a measurement 



of periastron advance for equal and unequal mass non- 
spinning black hole binaries. For eccentric binaries, peri- 
astron advance will result in a characteristic modulation 
of the observed GW signal, and hence it is important 
to quantify its frequency. We find that the numerically 
computed periastron advance fi$/51 r disagrees with both 
3PN predictions for equal-mass binaries, as well as with 
the test-mass limit of geodesic motion in a Schwarzschild 
background. As shown in Fig. [3 the periastron advance 
for black hole binaries lies roughly halfway between these 
two analytic calculations. The unequal mass evolutions 
considered have very small eccentricities; this is unfor- 
tunate for our current purposes, as this made it impos- 
sible to measure periastron advance well enough to test 
reliably the approach to the test-mass limit with increas- 
ing mass-ratio. While the data appear to approach the 
test-mass limit as the mass ratio deviates from unity, 
cf. Fig. [51 detailed confirmation will have to await until 
this analysis is repeated with somewhat higher eccentric- 
ity runs in the future. Nevertheless, even the equal-mass 
case shows that periastron advance is yet another feature 
of fully numerical calculations that is not accurately pre- 
dicted by post-Newtonian expansions. To achieve agree- 
ment, one may have to go to higher order post-Newtonian 
expansions, or one may have to incorporate finite-size 
effects. More pragmatically, for applications to grav- 
itational wave data-analysis, one might also introduce 
fitting parameters into the post-Newtonian models, and 
choose these parameters to enhance agreement with the 
numerical waveforms. 
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Appendix A: PN periastron advance 

In post-Newtonian approximations, theperiastron ad- 
vance was calculated to 3PN order in j47| for circular 
orbits in terms of the frequency-related parameter x. In 
the nonspinning circular case, the explicit expression for 
K is given by Eq.(5.11) of Ref. [47| in terms of the an- 
gular momentum density j for circular orbits and the 
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symmetric mass ratio v = miiri2 / (mi + TO2) 2 , where mi 
and m 2 are the masses of the two bodies, as 



1(45-12^)1 + 6 



rl35 Al 2 101. 



53 
— i 
24 



^static — ^ ^kinetic 



1 



(Al) 



where the value of the ambiguity parameter w s t a tic was 
computed by Ref. [52[ to be zero, and the ambiguity pa- 
rameter ^kinetic was shown to be 41/24 by Ref. [53|. The 
ratio 1/j 2 is replaced for circular orbits by l/j 2 irc where 



At the perihelia and aphelia of a test particle bound in 
an orbit around a black hole of mass M, r reaches its 
minimum r_ and maximum r+ when dr/dQ vanishes, so 
we can write 



1 A(r ± ) 
rl J 2 



E 

I 2 



(B5) 



From the above relation, the constants of motion E and 
J can be written as 



1 



xil 
16 



1 

" 3 
1 

64 



25 

(9 + vjx H vx 



41tt 2 - 



4 

5269 
~6~ 



511 
192 1 



3 / 2 
432^ _ (^static + ^ ^kind 



)]x 3 }.(A2) 



and 



E = 



A(r+)rl - A(r-)r 2 _ 



t2 A(r + )-A(r.) 



l/r\ 1/r 



(B6) 



(B7) 



Appendix B: Test-mass periastron advance for a 
Schwarzschild black hole 



Test particles follow geodesies in the background space- 
time, which here is given by the Schwarzschild metric: 



ds 2 



-A^^dt 2 + A(r)dr 2 +r 2 dn 2 , 



(Bl) 



where A(r) = (1 - 2M/r 



From Ref. j54[, these 



geodesic equations are given in term the radius r of the 
position vector as a function of time t by 



(B2) 



By integrating Eq. (|B4|) . we hnd that the angle swept 
out by the position vector as r increases from r_ to r+ 
is given by 



<f>(r + ) = <&(r_) 



'•+ 



A l l 2 (r) 



(B8) 

A(r) E 1 1 _1/2 dr 

J2 J2 r 2 r 2 



Then the orbit precesses in each revolution by an angle 
A$ defined as 



and 



A 2 (r) 



dr\ 2 J 2 
dt) + r 2 



A(r) 



(B3) 



where E and J are constants of motion. 

Since we are interested in measuring the periastron 
advance, we obtain the shape of the orbit using Eqs. (|B2[) 
and (lB3l): 



A(r) ( drX 
d&) 



1 



A(r) 
J 2 



E 
J 2 



(B4) 



A$ = 2|$(r + ) - $(r_)| - 2tt . 



(B9) 



To compute the periastron advance K as a function 
of the orbital frequency f}$, we pick a set of values for 
(r_ , r+ ) such that r + = r_ + e where e is a small positive 
number. The fractional periastron advance is estimated 
using Eq. (|B8|) . and the orbital frequency is estimated 
using Eq. (|B2|) . While if can be computed using elliptic 
integrals for Eq. (|B8|) . in practice it is simpler to evaluate 
it by numerical quadrature. 
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